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Introduction 



Motivations and a full account of Numerical Stochastic Perturbation Theory (NSPT) method 
for Lattice Gauge Theories, with special regard to dynamical fermions, can be found in ^ and 
references therein. In another contribution in this conference [|2|] we present our recent results. 
Here we want to stress some technical aspects, dealing in particular with fermionic observables 
computations. 



1. Langevin evolution 

According to the Stochastic Quantization approach, we randomly sample the phase-space of 
the generic field theory with action S[(p], according to the Langevin equation: 

t can be regarded as a new, non-physical, stochastic time. This sampling is such that the average 
over the gaussian noise X] leads, for infinite time evolution, to the same expectation values of 
Feynman's path-integrals: 

Perturbation theory is performed by thinking of the field as a formal power series in the coupling 
constant A: 

(j> n (x;t) — >£A>i n) (x;0. (1-3) 

n 

Plugging it in the Langevin equation results in a hierarchical system of differential equations, one 
for each perturbative term %' of the series. We solve it numerically via discretization of the 
stochastic time t = m. Every dynamical variable has to be replaced by its perturbative (truncated) 
expansion 

A ^ {AW}„<„ , (1.4) 

that is by a collection of dynamical variables. Consistently, every algebraic operation has to be per- 
formed 'order by order': for additions and scalar-multiplications this simply means to 'vectorize' 
them, 

(flA) w =flA(") (1.5) 
(A + B) (n) = AW + BW n = 0,l,...,n, (1.6) 

while ordinary multiplications have to be replaced by the (truncated) Cauchy product: 

(AB) W = £AWBW n = 0,l,...,n. (1.7) 

The need for a collection of variables in the place of a single one makes simulations very memory 
demanding. Moreover, since the computational effort of each basic operation is increased by the 
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order-by-order mechanism, the higher loop one wants to compute, the more expensive simulations 
are. In spite of this, two remarks are in order. First, there are cases where one can take advantage of 
the perturbative approach (see Sec]2|). The second remark is related to parallel-computing, where 
inter-nodes communications represent the bottleneck for a true scalability of performances (and 
hence of the overall computational power). The communications cost is proportional to the amount 
of data to be transferred, i.e. perturbative simulations communication time increases linearly with 
the maximum perturbative order n. On the contrary, floating point operations increase quadratically 
with n, namely as n(n — l)/2. So, as n increases, one improves the ratio between the computation 
and the communication times. 

2. Inverting the Dirac operator: unquenched dynamics and fermionic observables 

One of the main advantage of NSPT is that it is surprisingly cheap to invert the Dirac operator 
M. This is needed both for fermionic observables computations and unquenched dynamics. 
During dynamics, one has to compute the drift term of the Langevin equation 

V a Sp = -n f Tr [V a M M~ l ] (2.1) 

for the fermionic action = —n/Tr logM. The trick is to recover ( |Q| ) as an average over a 
gaussian noise t, 

\Sp = -n f [ $ (V a M) kI (M- 1 ) lj ^ ] , (2.2) 

so that one is left with the problem of inverting M over a source: A^.^Fj = In a perturbative 
context, higher orders *P can be recursively computed starting from the inversion of the 0-order 
only. The key point is that such a 0-order inversion is the tree-level (lattice) Feynmann propagator, 
which does not depend on the field configuration and is diagonal in momentum space. Hence, all 
the non-locality character of the problem is reduced to an iterative application (back and forth) of 
a (fast) Fourier transform (for details, see [[I]]). 

The same procedure is needed for the measurement of the quark propagator. Now we do not 
have to plug a noise £ but a 5-source for each element i of the propagator^ we are interested 
in: M k }¥j = 8[. A slight generalization of the same procedure can be used to compute fermion 
currents. The generic bilinear operator reads: 

8{M- l GM- l )8, (2.3) 

where G is a generic Dirac matrix: G = 1L, 75, y M , 7^75 and cfy v x [Yfi,7v]- It can be obtained by 
means of suitable scalar products involving again the inverse of Dirac operator: 

(8M- l )G{M- l 8) . (2.4) 

This strategy aims at saving CPU-time, but, as a drawback, it requires a lot of memory. Even if you 
are interested in (a small amount of) diagonal entry (in momentum space) of the current operators, 
you have to sum over the inner dummy indices, which, as said, extend all over the lattice volume. (2) 

Oft is worthwhile to note that we measure propagator in momentum space, as opposed to the non-perturbative 
practice. So in the 5-source, i is actually a momentum space index. 

®As a matter of fact, our APEmille machine does not have enough memory to perform such a computation on-the- 
fly, and so it has to break up the mechanism (carrying out some computations twice). On the other hand, instead, the 
memory equipment of a PC-cluster can do the job. 
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3. Exploiting the hypercubic symmetry to guide the continuum limit 

In order to fix ideas, let us focus on the two points vertex function T 2 . From a general point of 
view, with respect to Dirac-spin space, it can be decomposed into its component along the identity 
and the gamma matrices: 

r 2 = r id + ty = r id u + y r M . (3. i) 

n 

With respect to space-time, each component and can be (function of) Lorentz-Euclidean 
invariant quantities only. In the continuum 0(4) -symmetry case p 2 is the only scalar quantity at 
one's disposal and the only vector quantity. So ( |3.1[ ) would translate in 

r id ~ r id (p 2 ,m) (3.2a) 
r^ p ^T(p 2 ,m) i.e. y~/T(p 2 ,m). (3.2b) 

On the lattice one has to take into account the smaller hypercubic symmetry W4, the crystallo- 
graphic group of the discrete Tt/2 rotations on the six lattice planes of a 4-d hypercubic lattice onto 
itself, with the addition of the reflections. Less symmetry means less constraints, and so more as- 
sortment. Within the scalar sector, one has to face combinations of even powers of the momentum 
components: 

§ (2) Cp)=L/^ (3.3a) 

^\p) = M , ^\p) =11 py v , ^\p) =LLL pIpIpI ( 3 - 3c ) 

and so on. (3) Scalar invariants would suffice if we are interested in 1^ only, since, as said, itself is 
a scalar. , in turns, has to transform as a vector under the W4 symmetry, so that one has to take 
into account any odd (4) powers of p^ : 

^\P)=P, , ^\p)=PI , rfXp)=Pl , ••• (3-4) 
In the end the expressions (3.2 )'s for the Y 2 translate on the lattice into:^ 5) 



r Sd ~r id (sf ,) ) (3.5a) 
r M ~r M (sf } ,vi 2 " +1) ). (3.5b) 

By restoring physical dimensions, each lattice momentum gets its lattice-spacing factor a in 
front: p = ap. So one can recover the desired continuum limit a — > just in the limit of vanishing 



Strictly speaking, the ones listed above represent only a particular base for the generic invariants of a given order. 

2 



The equally invariant object (£„k] , for example, can be obtained by means of a suitable combination of the (3.3b )'s. 

Even powers are forbidden here — just likewise we considered, before, no odd powers for the § invariants — 
because the hypercubic group symmetry does include axes reflections. 

(5 'From now on we will consider the vanishing bare mass case, m = 0, so that all the mass-dependences can be 
neglected. 
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Figure 1: Lattice corrections and continuum extrapolation: example of scalar quantity (on the left) and 
'vector-like' quantity (on the right). On the left, interpolation fit is represented with circle around numerical 
data (with errorbars, almost invisible in the picture); on the right, the (more involved) interpolation fit results 
in the solid lines. See text for details. 



(lattice) momentum p. Therefore, a clear procedure for disentangling each quantity from finite 
lattice-size artifacts relies on a small-momenta expansion, so that only a few of them — the lower 



power ones — have to be really taken into account. As a result, the expansion for the (|3.5|)'s reads: 
T id = ff (0) + o^p 2 + V a,. (4) S (4) + V a^sf +..., (3.6a) 
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(3.6b) 



From dimensional analysis, the first coefficient of the series represents just the continuum limit 
value. Figure [l] shows two typical examples. The plot on the left represents measurements of (the 
third perturbative loop of) the critical mass for Wilson fermions. Since it is a scalar quantity, its 
behaviour is affected only by scalar lattice corrections, and so the picture is quite simpler with 
respect to the more general case (see next picture). For practical reasons data are shown versus 
the leading correction only, the momentum square, but the presence of the other invariants shows 
up in the jaggedness of the curve (which indeed does not represent statistical fluctuations of the 
signal). The plot on the right represents measurements of (the first perturbative loop of) the field 
renormalization constant, which results from the component along of r 9 . The situation is quite 
more involved with respect to the case of a scalar quantity, since here 'vector-like' lattice correc- 
tions mix up with the scalar ones. So we have not only the jaggedness of the scalar corrections, but 
points arrange into different levels according to the value of the component of the momentum in the 
direction along which the gamma-projection was performed. In simpler words: if the considered 
component of the momentum is equal to just one lattice-momentum-unit, then the associated point 
will fall on the lower level; if the component is twice, then the associated point will lie on the level 
above; if the component is three times the lattice-momentum-unit, then the point will lie on the 
next level; and so on. 
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4. Renormalization factors and anomalous dimensions 

In order to obtain the critical mass and the quark wave-function renormalization constant, one 
has only to get the two points vertex function inverting the fermion propagator. For bilinears, then, 
one uses the quark propagator (just measured for the same momenta) to amputate the external legs 
of the current operators. For UV-diverging quantities, care has to be taken while handling the 
logarithm terms, relying on the (already known) anomalous dimensions. 

As an example, let us refer to the scalar current in Landau gauge. The RF renormalization 
condition reads: 

Z q [ Z s (O s ) =1 (4.1) 

where Z q and Z s are the renormalization factors for the quark wave-function and the scalar current 
respectively, while O s is a convenient projection of the scalar current operator as it is measured on 
the lattice. Expanding everything in the lattice coupling j3 _1 up to, for instance, the first order, one 
gets: 

(0 (1) / 2^ (1) 



where is the one-loop anomalous dimension of the scalar current, which controls the ultra- 
violet divergence/ 6 ^ So, at one loop the renormalization condition becomes: 

-4 1) + ^ 1) -rf 1) logtf 2 ) + oP ) =0 (4-2) 



It is worthwhile to note the logarithm term of (4.2) just compensates the diverging behaviour of 



the lattice term o s , so that the whole expression is finite. Hence, the right prescription to get 
Zs^ is to subtract Y^hog(p 2 ) from Og , for each momentum, before taking the continuum limit. 
Actually, finite volume effects has to be taken into account as explained in [Q]. With these finite 
quantities, now, a little algebra (and the previous computed quark wave-function renormalization 
factors) allows to extract Zs . As order increases the algebra become somewhat more and more 
involved. In any case, at each order one is always able to extract the renormalization constant, in 
term of a finite number of log-subtractions. 
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^ Since we work in the Landau gauge, there is no anomalous dimension for the quark wave-function renormalization 
constant. 
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